home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
MATH
/
RMATH
/
CMATH.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1990-08-17
|
5KB
|
142 lines
{ unit CMATH/CMATH87; last modfied 17AUG90 }
{ Complex mathematics for Turbo Pascal }
{ Copyright 1990, by J. W. Rider }
unit {$ifopt N-} cmath {$else} cmath87 {$endif} ;
interface
uses {$ifopt N-} math {$else} math87 {$endif} ;
const maxcxray = {$ifopt N-} 5459 {$else} 3275 {$endif} ;
type complex = record rp, ip: float; end;
cxray = array [0..maxcxray] of complex;
function cabs (x: complex): float;
function phase(x: complex): float;
procedure conj (var x: complex); { complex conjugate }
procedure cneg (var x: complex); { x:=-x }
procedure cjx (var x: complex); { x:= j*x, j = sqrt(-1) }
procedure cnjx (var x: complex); { x:=-j*x, j = sqrt(-1) }
procedure crecip(var x: complex); { x:=1/x }
procedure csqr (var x: complex); { x:=x*x }
procedure csqrt (var x: complex); { x:=sqrt(x) }
procedure cln (var x: complex); { x:=ln(x) }
procedure cexp (var x: complex); { x:=exp(x) }
procedure csin (var x: complex); { x:=sin(x) }
procedure ccos (var x: complex); { x:=cos(x) }
procedure ctan (var x: complex); { x:=tan(x) }
procedure csinh (var x: complex); { x:=sinh(x) }
procedure ccosh (var x: complex); { x:=cosh(x) }
procedure ctanh (var x: complex); { x:=tanh(x) }
procedure casin (var x: complex); { x:=asin(x) }
procedure cacos (var x: complex); { x:=acos(x) }
procedure catan (var x: complex); { x:=atan(x) }
procedure casinh(var x: complex); { x:=asinh(x)}
procedure cacosh(var x: complex); { x:=acosh(x)}
procedure catanh(var x: complex); { x:=atanh(x)}
procedure cadd (var x,y: complex); { x:=x+y }
procedure csub (var x,y: complex); { x:=x-y }
procedure cmult(var x,y: complex); { x:=x*y }
procedure cdiv (var x,y: complex); { x:=x/y }
procedure cpow (var x,y: complex); { x:=pow(x,y)}
procedure cpoly (var x:complex; degree: integer; var coeffs);
implementation
{ "alphbeta" is used by the "cacos" and "casin" routines. I
derived it from Abramowitz and Stegun, Handbook of Mathematical
Functions. (as well as most of the routines) }
procedure alphbeta(var x:complex); var z:xfloat;
begin with x do begin
z:=system.sqrt(sqr(rp+1)+sqr(ip));
rp:=system.sqrt(sqr(rp-1)+sqr(ip));
ip:=(z+rp)/2; rp:=ip-rp; ip:=ln(ip+sqrt(sqr(ip)-1)) end end;
function cabs; begin cabs:=hypot(x.rp,x.ip) end;
procedure cacos; begin alphbeta(x); x.rp:=acos(x.rp) end;
procedure cacosh; begin cacos(x); cjx(x) end;
procedure cadd; begin x.rp:=x.rp+y.rp; x.ip:=x.ip+y.ip end;
procedure casin; begin alphbeta(x); x.rp:=asin(x.rp) end;
procedure casinh; begin cjx(x); casin(x); cnjx(x) end;
procedure catan; var rp2,ip2,z:xfloat; begin
rp2:=sqr(x.rp); ip2:=sqr(x.ip); z:=atan(2*x.rp/(1-rp2-rp2))/2;
x.ip:=ln((rp2+ip2+2*x.ip+1)/(rp2+ip2-2*x.ip+1))/4;
x.rp:=z end;
procedure catanh; begin cjx(x); catan(x); cnjx(x) end;
procedure ccos; var z:xfloat; begin
z:=cos(x.rp)*cosh(x.ip);
x.ip:=-sin(x.rp)*sinh(x.ip); x.rp:=z end;
procedure ccosh; var z:xfloat; begin
z:=cosh(x.rp)*cos(x.ip);
x.ip:=sinh(x.rp)*sin(x.ip); x.rp:=z end;
procedure cdiv; var w,z:xfloat; begin
w:=sqr(y.rp)+sqr(y.ip); z:=(x.rp*y.rp+x.ip*y.ip)/w;
x.ip:=(x.ip*y.rp-x.rp*y.ip)/w; x.rp:=z end;
procedure cexp; var z:xfloat; begin
z:=exp(x.rp); x.rp:=z*cos(x.ip); x.ip:=z*sin(x.ip) end;
procedure cjx; var z:xfloat; begin z:=-x.ip; x.ip:=x.rp; x.rp:=z end;
procedure cln; var z:xfloat; begin
z:=ln(cabs(x)); x.ip:=phase(x); x.rp:=z end;
procedure cmult; var z:xfloat; begin
z:=x.rp*y.rp-x.ip*y.ip; x.ip:=x.rp*y.ip+x.ip*y.rp; x.rp:=z end;
procedure cneg; begin x.rp:=-x.rp; x.ip:=-x.ip end;
procedure cnjx; var z:xfloat; begin z:=x.ip; x.ip:=-x.rp; x.rp:=z end;
procedure conj; begin x.ip:=-x.ip end;
procedure cpoly;
{ coeffs is assumed to be an array [0..degree] of "complex" }
var y:complex;
begin y.rp:=0; y.ip:=0; while degree>=0 do begin
cmult(y,x); cadd(y,cxray(coeffs)[degree]); dec(degree) end;
x:=y end;
procedure cpow; begin cln(x); cmult(x,y); cexp(x) end;
procedure crecip; var w:xfloat; begin
w:=sqr(x.rp)+sqr(x.ip); x.rp:=x.rp/w; x.ip:=-x.ip/w end;
procedure csin; var z:xfloat; begin
z:=sin(x.rp)*cosh(x.ip); x.ip:=cos(x.rp)*sinh(x.ip); x.rp:=z end;
procedure csinh; var z:xfloat; begin
z:=sinh(x.rp)*cos(x.ip); x.ip:=cosh(x.rp)*sin(x.ip); x.rp:=z end;
procedure csqr; var z:xfloat; begin
z:=sqr(x.rp)-sqr(x.ip); x.ip:=2*x.rp*x.ip; x.rp:=z end;
procedure csqrt; var absx,z:xfloat; begin
absx:=cabs(x); z:=sqrt((absx+x.rp)/2);
x.ip:=sgn(x.ip)*sqrt((absx-x.rp)/2); x.rp:=z end;
procedure csub; begin x.rp:=x.rp-y.rp; x.ip:=x.ip-y.ip end;
procedure ctan; var z:xfloat; begin
z:=cos(2*x.rp)+cosh(2*x.ip);
x.rp:=sin(2*x.rp)/z; x.ip:=sinh(2*x.rp)/z end;
procedure ctanh; var z:xfloat; begin with x do begin
z:=cosh(2*rp)+cos(2*ip); rp:=sinh(2*rp)/z; ip:=sin(2*ip)/z
end end;
function phase; var z:xfloat; begin z:=atan2(x.rp,x.ip) end;
end.